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Abstract The number of publications of aperture- synthesis images based on optical long- 
baseline interferometry measurements has recently increased due to easier access to visi- 
ble and infrared interferometers. The interferometry technique has now reached a technical 
maturity level that opens new avenues for numerous astrophysical topics requiring milli- 
arcsecond model-independent imaging. In writing this paper our motivation was twofold: 1) 
review and publicize emblematic excerpts of the impressive corpus accumulated in the field 
of optical interferometry image reconstruction; 2) discuss future prospects for this technique 
by selecting four representative astrophysical science cases in order to review the potential 
benefits of using optical long baseline interferometers. 

For this second goal we have simulated interferometric data from those selected astro- 
physical environments and used state-of-the-art codes to provide the reconstructed images 
that are reachable with current or soon-to-be facilities. The image reconstruction process 
was "blind" in the sense that reconstructors had no knowledge of the input brightness distri- 
butions. We discuss the impact of optical interferometry in those four astrophysical fields. 
We show that image reconstruction software successfully provides accurate morphological 
information on a variety of astrophysical topics and review the current strengths and weak- 
nesses of such reconstructions. 
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We investigate how to improve image reconstruction and the quality of the image pos- 
sibly by upgrading the current facilities. We finally argue that optical interferometers and 
their corresponding instrumentation, existing or to come, with 6 to 10 telescopes, should be 
well suited to provide images of complex sceneries. 

Keywords Instrumentation: interferometers • Techniques: interferometric • Methods: data 
processing • planetary systems: protoplanetary disks • stars: supergiants • stars: AGB - 
circumstellar matter • stars: imaging • Galaxies: active • Galaxies nuclei 



1 Introduction 



The recent years have seen a significant increase in scientific publications making use of 
images based on homodyne optical long-baseline interferometr y (e.g. |Monnier et al[|2007| 
Le Bouquin etal| [20091 |Schmitt et al[ [20091 |Lacour etaIl|2QQ9l|Kraus et al[|2009||Haubois| 



et al||2009| |Kloppenbor g"et al||2010[ |Millour et al| |2011| ). While model fitting of visibilities 



remains today the most accurate way to characterize an object brightness distribution, the 
technique has now reached a technical maturity level that opens new avenues for numerous 
astrophysical topics requiring milli-arcsecond model-independent imaging. 

These results are the culmination of years of efforts that started with speckle interfer- 
ometry, aperture masking and long-baseline interferometry. For the latter active or defunct 
prototype facilities such as COAST, IOTA, I2T, GI2T, Mark III, NPOI, PTI, SUSQ have 
progressively paved the way for large diameter telescope arrays such as CHAR A, Keck 
Interferometer, and VLTI. 

In order to assess the main astrophysical questions that visible and infrared interfer- 
ometry can actually answer, we first review the status of optical interferometry in terms of 
image synthesis achievements and then follow up with a discussion on the important issues 
from four main astrophysical fields namely planetary formation, stellar atmospheres, cir- 
cumstellar shells, and active galactic nucleQ Then, we simulate interferometric data from 
those environments and provide the reconstructed images that are achievable with current or 
near-future facilities, in order to discuss the impact of such techniques in those fields. 

The paper is organized as follows. The principle of interferometric imaging as well as 
a brief description of present imaging interferometers is given in Sect. [2] The four astro- 
physical topics are briefly described in Section [3] and for each one the parameters of the 
simulations are presented. The imaging facilities are reviewed in Section [4] In Sect. [5] we 
describe the process which allows us to produce realistic interferometric data and then re- 
constructed images. Finally, Section [6]is devoted to a discussion of the results and the future 
perspectives for optical long-baseline interferometry. 



2 Imaging with a long baseline interferometer 

2. 1 Principles of optical long baseline interferometry 

Principles of long baseline interferometry are well explained in numerous articles, confer- 
ence proceedings and monographies. The reader is referred to Thom pson et allpOOl) for a 

1 ISI belongs to the same fraternity but operates in heterodyne mode and will therefore not be discussed. 

2 This choice, obviously not exhaustive, should not hide the wealth of astronomical topics requesting milli- 
arcsecond resolution imaging: Cepheids, magnetic, Be, O, supermassive stars, stellar mass loss, jet formation, 
dynamics of close stellar clusters, SMBH galaxies etc. 
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detailed description of the technique which lays out nicely most of the fundamental prin- 
ciples, albeit mostly oriented towards radio astronomy. More recently there have been in- 
teresting collective works, review articles and monographies dedicated to the specifics of 
optical long baseline interferometry (e.g. |Lawson| [2000} |Monrnerl [20031 [HanifT| |2007] |MaT| 
|bet a nd Perrin, 2007 ; Glindemann, 201 1}. In our definition, optical interferometry covers the 
wavelength regime from the visible 0.4 jim to the mid infrared 10-20 jim. 

An interferometric array measures the spatio-temporal coherence of the electromagnetic 
field through two or more telescopes (e.g. Goodman, 1985). The core observable of an in- 
terferometer is the complex visibility Vja' 



V a = \V a \a&*» (1) 

Where |Vy is the amplitude of the visibility and &u its phase. Vy is extracted from the 
measurement of the contrast and phase of the interferogram formed between two telescopes 
k and /. Unlike radio interferometers, where the detection of the electromagnetic field is 
done at the telescope level (heterodyne interferometry), most optical interferometers require 
the actual formation of interferometric fringes on a detector. 

The fundamental theorem of long baseline interferometry is the Van-Cittert Zernike the- 
orem ( [Good man , 1985 ; Thomp son et al[|200T] }. It relates the complex visibility V measured 
at a spatial frequency B/A (B being the telescope baseline vector and A the wavelength) to 
the actual object brightness I(x) distribution through a Fourier transform relation. 



f°° /(x)exp" 2 ^' xBM dx 
V(B/A) = J -°° f00 (2) 

The components (u, v) of B/A form the spatial frequency coordinate system, x represents 
the angular position vector. Therefore, interferometric observations provide information 
about spatial frequency components of the brightness distribution. Similarly to the "Rayleigh 
criterion" used in single pupil telescopes, one can define the resolution of a two-telescope 
interferometer of baseline B as A/2B (in radians) which corresponds approximately to the 
milli-arcsecond level in the near infrared with typical hectometric baselines. The amplitude 
of the visibility is related to the projected brightness angular size while the phase provides 
information on the brightness photocenter location. 

In practice, visible and infrared interferometry require mixing the light received from an 
astronomical source and collected by several independent telescopes separated from each 
other by tens or even hundreds of meter^] The light beams are then overlapped and form 
an interference pattern if the optical path difference between the different arms of the inter- 
ferometer — taking into account paths from the source up to the detector — is smaller than 
the coherence length of the incident wave (typically of the order of several microns). This 
interference pattern is composed of fringes, i.e. a succession of stripes of faint (destructive 
interferences) and bright (constructive interferences) intensity. By measuring (i) the contrast 
of these fringes, i.e. the normalized flux difference between the maximum and minimum 
intensity, also called the visibility amplitude, and (ii) their phase, i.e. their position, one can 
construct an estimation of the so-called complex visibility. 



With the notable exception of the heterodyne ISI interferometer ( Townes and Wishnow 2008 
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2.2 Obstacles to imaging with an optical interferometer 

In principle, measuring the complex visibilities at all spatial frequencies up to a maximum 
value uv max , should permit the retrieval of the synthesized image of the spatial intensity 
distribution of the object with a spatial resolution l/uv max by inverse Fourier transforming 
the 2-D visibility map. However there are two main obstacles to this simple inversion. 

First, in all practical cases, the (u, v) plane cannot be fully covered, because the spa- 
tial frequencies sampled by an interferometer are limited by the number of pairs of tele- 
scopes. With the Earth rotation these spatial frequencies follow (u, v) tracks in shape of 
arcs of ellipses (Thompson et al, 2001; Segransan] |2007| ) which provides super-synthesis 
(also called "Earth rotation synthesis") and therefore increases the coverage. In the partic- 
ular cases where no strong wavelength dependence of the object brightness distribution is 
expected, one can use measurements at different wavelengths to help further extend the (w, v) 
coverage, but only radially. 

Secondly, the phases of the visibilities are frequently lost because they have been scram- 
bled by atmospheric blurring effects ( Quirrenbach , 2000). Instrumental effects also con- 
tribute but in a more static way. In the optical domain (as in the radio but on a wilder 
scale), rapid fluctuations in the atmospheric optical index induce random optical path vari- 
ations. Typical atmospheric coherence times at these wavelengths are in the range of 1 to 
10 milliseconds. Consequently, measuring direct phase information is challenging. Several 
techniques have been explored to circumvent this difficulty, all of which have inherited from 
experience in radio interferometry, and are well detailed in Monnier (2007|). 

- The phase of the product of the 3 complex visibilities measured with 3 telescopes form- 
ing a closed triangle, called the closure phase, was first introduced by Jennison ( 1958 ) 
in order to cancel out instrumental and atmosphere-induced phase errors and recover 
partial phase information. It was later pioneered in the optical domain by Baldwin et al 
( 1986 ) and subsequent aperture masking experiments at Cambridge and is now a routine 
observable of optical interferometers. Parallel to that, triple-correlation techniques used 
in speckle interferometry (e.g Weigelt, 1991 ) also make use of closure phase although 
the formal relation between the two techniques was established by Roddier (1986). In 
practice, the closure phase is estimated from the phase of the so-called "bi-spectrum" 
Bijk which is constructed from the product of complex visibilities estimated in a closed 
triangle B ijk = Vy V jk V ki . 

- The self-calibration method as described by |Cornwell and W ilkinson] ( |1981| ) aims at 
retrieving phase information on a maximum number of baselines from closure phase 
estimators. It integrates the phase recovery procedure in the image reconstruction pro- 
cess. For that purpose it uses an iterative scheme that starts from a trial image and seeks, 
through modeling and fitting of atmospheric and telescope phase errors, to find the best 
phase estimation consistent with closure phases. 

- Phase referencing encompasses various techniques that aim at measuring, as directly 
as possible, the phases. The idea is to calibrate the phase measurement with a phase 
reference, either internal or external. Differential phases, which are obtained by com- 
paring phases at different wavelengths with a reference channel, provide direct means 
of phase estimation. This is true provided one can constrain the morphology at the ref- 
erence wavelength (e.g an unresolved emission in a continuum channel) and model the 
atmospheric and instrumental differential phase. The reader is referred to the work of 
Vakili et aT|( 1997|), for example, for illustration of the differential phase extraction, and 



Schmitt et al ^2009 ) and previous NPOI team publications for a detailed demonstration 
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of spectral differential imaging. Astrometric phase referencing uses the capability of ob- 
serving two objects within the atmospheric isoplanatic patch to provide an astrometric 
reference to the phase (with the motivation of authorizing longer integration times if a 
nearby bright source can be used as a reference). This requires a specific infrastructure 
allowing for dual-star observation. It has been implemented at PTI, NPOI, and is being 
implemented at Keck and VLTI (e.gjpelplancke, 2008 ) but has not led yet to notable 
mapping results. 

2.3 Principles of image reconstruction 

Because of the voids in the (u, v) plane coverage and the lack of complete phase information, 
a given set of data can be fitted by several different brightness distributions. The image power 
spectra and bispectra of all these distributions fit the data within the error bars in the least 
squared sense. Hence, additional prior information is required to discriminate between those 
images and to further constrain the solutions. The best image is defined as the most probable 
image given the data and prior information. The bias toward the prior is the price to pay to 
select an image from such sparse data. Reconstructing an image consists then in minimizing 
a quantity connected to the data and to the prior with respect to the pixel values of the image 
jThiebaut and Giovanneim|2010] >: 

Xbest = argmin x€iQ [/ da ta(*; j) + V /priori)] (3) 

where x are the image parameters (e.g. the pixel values) constrained to belong to the feasible 
set Q of, e.g., non-negative and normalized images, /data measures the discrepancy between 
the image model and the data y, f pr \ or is the regularization term which enforces the priors 
and the hyperparameter ji is a weight factor that tunes the balance between these two terms. 

The term /data is usually derived from an analytical model of the data and the statistics of 
the noise, although slightly different approximations (e.g., Meimon et al , 2005 ) are made by 
current image synthesis algorithms. Note that /data m ay be itself the sum of different terms 
to account for different kinds of data such as power spectra and phase closures. 

The term / pr j 0r is used to inject in the minimization process a priori knowledge on the 
object. In addition to fundamental constraints such as positivity (i.e. intensity can not be 
negative) and limited field of view there are numerous possible priors that can be included 
in the regularization. Although automation of the image reconstruction process is highly 
desirable, most of current (u, v) coverages in optical interferometry leave room for multiple 
possible solutions. Therefore the observer cannot avoid introducing, through the regulariza- 
tion, his knowledge of the expected brightness distribution. The interested reader is referred 
to the recent work by Renard et~al| ( |201 1| ) for a benchmarking of regularizations. 

There are several popular regularizations. 

1 . Maximum Entropy Method (MEM, |Gull and SkilEngl [T984| |Skilling and Bryan| |1984a[ 
Nara yan and N ityananda] |1986| ). This is one of the most employed techniques. It tries 
to find the smoothest image compatible with the data while keeping high frequency 
information. For that purpose it uses a smoothness scalar estimator called entropy which 
has several possible definitions, one of the most popular being the Gull-Skilling entropy 
(Skilling and Bryan, 1984b). MEM will find the image that best fits the data and keeps 
entropy maximum. Because of the smoothness requirements MEM will introduce super- 
resolution to the image, i.e. it will put information at spatial frequencies higher than the 
usual optical Rayleigh criterion. 
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2. Compactness allows the definition of a spatial zone outside of which no flux is expected. 

3. Total Variation tries to minimize the total gradient of the image. It favors parts of the 
image with uniform zones and fast intensity changes (like a stellar surface). It has been 
found by Renard et al (2011) as one the best "all-terrain" regularization methods when 
applied to different object morphologies. 

4. Smoothness regularization smooths the image but can be adjusted to preserve some sharp 
variations at edges or point-like structures, favor some regions in the image etc. There 
are various ways to force and control smoothness. 

5. Spars eness aims at restoring an image with the smallest number of components {e.g. 
CLEAN method in radio-astronomy or Building Block Method in optical interferometry, 
see below). 

Figure [TJillustrates visually the effect of different prior regularizations on different types 
of brightness distributions. The first column is a doughnut like structure surrounding a cen- 
tral star. The middle column is a photosphere with limb-darkening. The last column is a 
galaxy composed of extended and point-like emission. The prior regularizations, from top 
to bottom, are: the original image, smoothness, compactness, total variation, smoothness 
with sharp boundaries (e.g I2 ~ A norm), l p norm with p - 0.5 and p - 2, squared-root 
MEM and MEM with gaussian a priori image. 



2.4 Imaging algorithms and their application 



In the radio domain the CLEAN algorithm (Hog bom||1974 » and its different flavors asso- 
ciated or not with self-calibration ( [Corn well and W ilkinson, 1981] ) have been the dominant 
tools for image reconstruction. This is not the case in optical long baseline interferome- 
try where mapping techniques are still the subject of active research. It is essentially the 
absence of direct phase information that has triggered new developments with respect to 
radio-astronomy. Although this effort started with the rise of aperture masking (Baldwin 
|et al[|198 6; Haniff et al, 1987; Readhead et al, 1988) in single-dish telescopes, the similari- 
ties are so strong that both techniques can use the same tools. For our purposes we can safely 
ignore the specificities. Here we discuss the methods that are being actively used in actual 
astrophysical programs and are sometimes tailored to deal with different source structures 
and the number of telescopes used. 



The Building Block Method (BBM) has been proposed by Hofmann and Weigelt ( 1993} and 
offers many similarities with the CLEAN method except that it can handle bispectra (and 
hence closure phases). The minimization is carried through a matching pursuit algorithm 
that imposes the sparsest solution. 



The BiSpectrum Maximum Entropy Method (BSMEM) uses the MEM method to regular- 
ize the image reconstruction process using visibility amplitudes and bispectra. It supposes 
independent gaussian statistics for amplitude and closure phase and uses the Gull-Skilling 
( [Skilling and Bryan) |1984b| ) entropy as regularization criterion. Image entropy provides a 
metric to measure smoothness in the image. This criterion will minimize the regularization 
weight based on the comparison with a reference image that can be updated through an iter- 
ative process. It uses a specific numerical optimizer MEMSYS that implements the strategy 
established by Skilling and Bryan ( 1984a) and can handle any type of data. 
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Fig. 1 Effects of different priors. From left to right: reconstructed images of a doughnut like structure sur- 
rounding a central star, a photosphere with limb-darkening and a galaxy composed of extended and point-like 
emission. From top to bottom: original image, smoothness, compactness, total variation, smoothness with 
sharp boundaries (e.g norm), l p norm with p = 0.5 and p = 2, squared-root MEM and MEM with 

gaussian a priori image. 
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Table 1 Comparison of the image synthesis algorithms used in optical interferometry. 



Code 


Data 


Regularization 


Optimization Strategy 


name 








MiRA 




pUblllvliy, lUld.1 VdlltlllUIl, €2? 

h-h, Gull-Skilling entropy 
with floating or given prior 
image 


111111 ICQ lllClllUiy CjUdal 

Newton with bound (pos- 
itivity) and normalization 
constraints 


BSMEM 


bispectra, power spectra, 
complex visibilities 


Gull-Skilling entropy with 
given prior image, multi- 
scale entropy 


non-linear conjugate gradi- 
ents with unsupervised hy- 
per parameter control 


WISARD 


pseudo-complex visibili- 
ties formed from phase 
closures and power spec- 
tra 


positivity, £2, 


alternative minimization and 
self calibration 


BBM 


bispectra 


positivity, sparseness 


matching pursuit 


MACIM 


any 


any 


global optimization by sim- 
ulated annealing 



MEM and self-calibration Here, the maximum entropy method is embedded in an interac- 
tive process that starts with a guess on the instrument phases compliant with the closure 
phase measurements. The iteration includes fitting telescope atmospheric and instrumental 
phase errors and remapping of the image while respecting the closure phases. It therefore 
iterates in order to converge to the most probable image using the best guessed phases. 



The MACIM algorithm (MArkov Chain Imager Ireland et al, 2006), unlike the preceding 
tries to find globally the images that are optimal in the Bayesian sense (and not the image 
that apparently fits best the data). For that purpose it tries to maximize the a posteriori 
probability: 



Pr(x, y) = exp 



2 (/data ' 



" M/prior) 



(4) 



and therefore can provide a joint probability density of the probable images. It can use or 
not use regularizations (currently MEM and dark zone connectivity) and handles any kind 
of data. 



The MIRA algorithm (Multi-aperture Image Reconstruction Algorithm, Thiebaut, 2008 ) 
minimizes Eq. ([5]) by using a non linear optimization algorithm. Because the method does 
not allow a global optimization, it will depend on using the initial image as a constraint. 
It does not try to reconstruct the phases directly and therefore can handle closure phases 
data and any other observable. It encodes many types of regularization which offers the 
possibility of testing the reliability of the reconstructed images. 

The WISARD algorithm (Weak-phase Interferometric Sample Alternating Reconstruction 
Device, Meimon et al, 2008), is an application of the self-calibration idea to optical inter- 
ferometry. It reconstructs the phases from the closure phases. The phase inversion process 
has multiple minima and therefore the global minimization is a multi-modal problem and 
can depend on the initial image. 

Table [T] summarizes the features of the different algorithms that were specifically de- 
signed to cope with optical interferometry data. These algorithms are the main challengers 
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of the Interferometric Imaging Beauty Contests which, every two years since 2004, quanti- 
tatively compares the results of various image synthesis methods on simulated optical inter- 
ferometric data ( [Lawson e t al, 2004 1 [2006} |Cotton et al] |2008l |Malbet et al|[20To|). Among 
these algorithms, BSMEM ( |Buscherl|1994[lBaron and Ymmg||20081 > and Mi^A JThiebautl 
2008} have repeatedly won the contest for data sets typical of next generation instruments. 
These are the ones that were used to produce the images from simulated data sets presented 
in this paper, together with MACIM ( [Ireland etal||2006| ). 



2.5 Examples 

In this section, we describe some selected significant astrophysical publications making use 
of reconstructed imaging (Fig. [5}. While we acknowledge that these choices suffer from 
selection effects, (for example we have limited the amount of binary images) we believe that 
these examples illustrate well the state of the art. 

Image 1 in Fig.[2]shows one of the early attempts by Quirrenbach et al (1994 ) to recon- 
struct images out of visibility amplitude data using the Mark III interferometer. Even with- 
out phases they were able, using MEM-inspired methods, to map the H a emission around £ 
Tauri showing clearly its extension and flattening (related to its inclination). The absence of 
phases appears in the centro- symmetry of the image. 

Images 2 and 3 in Fig. [2] show the pioneering images of binaries obtained at COAST 
( [Baldwin ef al, 1996) and NPOI (now NOI) ( [Benson etal[|1997| ) using standard radio pack- 
ages associated with self-calibration (for the NOI reconstruction). 

Young et al (2000 ) used aperture masking and the COAST to map the surface of Betel- 
geuse. In image 4, the aperture masking image obtained with the William Herschel Tele- 
scope is shown. The surface is marginally resolved since the angular resolution is barely 
sufficient. The image was obtained with the combination of the MEM method and self- 
calibration. More recently Haubois et al (2009 ) imaged the spotty surface of Betelgeuse 
using IOTA (images 5 and 6 of Fig. [2}. Both WISARD and MIRA reconstructions lead to 
the same result confirming what had already revealed by Young et al (2000 ). In this case, 
both the dynamic range and the ability to characterize the spots are limited by the instrument 
accuracy and the difficulty for an interferometer to lock on low visibilities. This constraint 
has to be taken into account in the design of any future "imaging instrument", as we will 
discuss later. 

Images 7 to 9 of Fig. [2] are reconstructed images of respectively Altair, T Leporis and 
Epsilon Aurigae. Images 7 and 9 were obtained with CHARA using MACIM and BSMEM 
(respectively Monnier et al, 2007; Kloppenb org et al[[2010| ). Image 8 was taken with VLTI 
using MIRA ( |Le Bouquin et al[ 2009 \ . These images are an illustration of how image recon- 



struction can retrieve moderately complex morphologies with a limited (u, v) coverage that 
might have escaped traditional model fitting strategies. 

Spectral phase referencing has proven to be a powerful tool to image line emitting re- 
gions. For this purpose, a sufficient spectral resolution is required in order to compare purely 
continuum emission to line emission interferometric observables. Results using this tech- 
nique are illustrated in images 10 and 11 of Fig. [2] where the H a emission around the star 
ft Lyrae has been mapped using the AIPS package ( Sch mitt et al[ |2009| ). This was possi- 
ble because each baseline phase corresponding to the line emission was isolated from the 
continuum. In a similar spirit, Mill our et al| ( |2011| ) have used AMBER with high spectral 
resolution, MIRA, and a self-calibration approach to retrieve absolute phases in the Br r line 
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of the supergiant HD 62623. Since the line was both spectrally and spatially resolved it was 
possible to reveal the kinematics of a disk around the central star. 

Mapping the infrared emission of protoplanetary disks has always been one of the key 
science drivers for optical arrays and their focal instrumentation. Images 12 to 14 of Fig. [2] 
are an illustration of the best disks maps achieved so far in this field. Images 12,13,14 have 
been obtained, respectively, by Aperture Masking at Keck (jTuthill et~al] |2002| ) and long 
baseline interferometry at VLTI ( |Kraus et al[ [201^ |BenTsty et al[|20fTT Tlmages 12 and 14 
clearly show the inner rim of the circumstellar disk caused by dust sublimation. Maps 13 
(using the BBM method) and 14 (using MIRA) display limited features in the image because 
of the relatively sparse (u, v) coverage obtained with the VLTI. 

The images obtained with aperture masking at Keck by Tuthill et al ( 1999, Wolf-Rayet 
104) and Monnier et al (2004, Evolved star NML Cyg) show the most complex features of 
all the images selected here (images 15-17 of Fig. [2]). This complexity is a direct illustration 
of the importance of a sufficient (u, v) coverage since both used masks with 21 apertures. 
MEM priors were chosen in the reconstruction process. The case of the NML Cyg map 
is also very illustrative of the requirement to strategically pave the (u, v) plane. Therefore, 
one should always consider complementing long-baseline interferometry with single-dish 
aperture masking if possible or vice-versa. In this particular case, the Keck aperture masking 
data lacked high spatial frequency constraints on the central source (left image). The addition 
of IOTA data allowed the central source extension to be characterized and improved the 
quality of the reconstruction. Conversely, extended fluxes might alter interferometric data 
and should therefore be at least calibrated and ideally measured with short spatial frequency 
instruments such as aperture masking instruments. 



3 The future of aperture synthesis: as trophy sical objects at the milli-arcsecond scale 

While optical long baseline interferometers have made tremendous technical and scientific 
progress, we wanted to evaluate what improvements could be brought to the field of aperture 
synthesis imaging. Following the recommendations of the SOC supervising the Workshop on 
Interferometric Imaging 2009 (see acknowledgements), we selected four main astrophy sical 
topics for which visible and infrared interferometers are expected to bring breakthroughs 
and for which milli-arcsecond scale imaging will certainly improve our understanding of 
the physical phenomena. We restricted ourselves to the study of the continuum emission of 
various object classes. The following sections detail those topics. 



3.1 Planetary signatures in protoplanetary disks 
3.1.1 Context and major issues 

Currentl}0 more than 700 exoplanet^] together with several thousand of planet candidates 
have been detected so far and although the detection techniques used may cause strong se- 
lection effects, it is now obvious that the solar-system planets cover a small region in the pa- 
rameter space of possible planetary characteristics. To understand the observed diversity of 
planets, one has to consider the history of the planetary systems by studying their formation 



4 as of December 201 1 

5 http://exoplanet.eu 
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Fig. 2 The wall: examples of published reconstructed images (see text for details). 1 Quirrenbach et al (1994 » 
2 Bal dwin et aT]fl 996 1 3 Be nson et al]dl997) 4|Young et alH2000) 5&6 Haubois et al (2009) 7 Monnie r et al 
(2007|,8: |Le Bouquin et al| 12009), 9[K loppenbor g~et al|(|2010),10|Scnmitt et alH2009>, 11 : |Millour et al 
j20TTj,12|Tuthill et alH2002) rT3lKraus et al (2010», 14 Benis ty et alH2011) , 15 |Tuthill et alH1999) , 16&17 
|Monnier et alH2004) . Credit for figure 8: ESO/J.-B. Le Bouquin et al. 
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and evolution. Triggered and supported by the rapidly growing field of the search for extra- 
solar planets and the observational material already collected, our theoretical understanding 
about the planet formation process has been improved significantly. Nevertheless, the link 
between theoretical models and observations is still rather weak. This link is provided by 
circumstellar disks in which planets are expected to form. These disks are considered as the 
natural outcome of the protostellar evolution, at le ast in the case of low and intermediate 
mass stars (e.g. | Adams et al 1987 Lissauer 1993|). These disks have masses (10~ 3 - 1 M ; 
Beckwith jTal||199Q| ) andsizes (50 - 1000 AU; |McCaughrean and O'DelllfT^ |Dutrey| 
et al |1996| ) comparable to those expected for the primitive solar nebula. 



Theoretical investigations show that the planet-disk interaction creates structures in 
young circumstellar disks which are usually much larger in size than the planet itself and 
thus more easily detectable (Masset and Papaloizou, 2003; Varniere et al, 2004; Crida and 
|Morbidelli] |2007| ). The specific result of the planet-disk interaction depends on the evo- 
lutionary stage of the disk. Typical signatures of planets embedded in disks are gaps and 
spiral density waves in the case of young, gas-rich protoplanetary disks, and characteristic 
asymmetric density patterns in older debris disks. Numerical simulations demonstrate that 
high-resolution imaging performed with observational facilities (already existing or soon 
becoming available) will allow these "fingerprints" of planets to be traced in protoplane- 
tary and debris disks (e.g. [Wolf and D'Angelo] |2005| ) but the prospect of such detections 
lies mostly in the mm domain. These observations will provide a deep insight into specific 
phases of the formation and early evolution of planets in circumstellar disks. For example, 
| Wolf and Klahr| ( |2002| ) showed that the Atacama Large Millimeter Array (ALMA) will al- 
low tracing of hyper-dense vortices which are supposed to represent an early stage of planet 
formation. Some authors have tested the detectability of planet-induced features, such as 
temperature increase due to planetary accretion, in the near- to-mid-infrared domain with 
interferometers (e.g. [Wolf, 2008) and have already pointed out that this represents a major 
challenge. These studies have to be pursued. 

For this work we wanted to test the limit of the planet-induced structures that could be 
detected. We have therefore not considered other crucial observational diagnostics that could 
also (and already do) allow the properties of protoplanetary environments at the astronomi- 
cal unit scale to be further constrained. These observable diagnostics include the disk inner 
rim vertical structure, the dust content and its vertical distribution, kinematics, asymmetries 
etc. Since these are already the subject of active research we instead chose to investigate 
more challenging observables which could become accessible in the future. 



3.1.2 Simulated images of disks perturbed by planets 



The model for radiative transfer in the vicinity of an embedded protoplanet used in this work 
is detailed in |Jang-Co ndell (2008 ). It is based on a one-dimensional analytical solution for 
radiative transfer in an optically thick protoplanetary disk ( [D'Alessio et al[[T998] ) extended 
to three dimensions. The calculation is carried out iteratively for self-consistency between 
temperature and density perturbations. The generation of simulated images is described in 
|Jang-Condelll ( [2009| ). 

The synthetic image produced for the image reconstruction is of a face-on disk at a 
distance of 140 pc. Since the simulation boundaries are limited in radius at both the inner and 
outer edges, the disk brightness is smoothly extrapolated as a power-law interior and exterior 
to the simulated region in order to eliminate image edge artifacts. The star is modeled as a 
black body of radius R* = 2.6 R Q and effective temperature T = 4280 K, corresponding to a 
1 Myr old star of mass M* = 1M ( Sie ss et al[|2000| ). 
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Fig. 3 Simulated images of our 4 scie nce cas es. Top, left to right: (a) disk in the N band with shadows created 
by low-mass embedded planets (Sect. |3.l!2) , inset shows disk with gap in the N band; (b) stellar surface of 



a supergiant star with convective cells in the K band (Sect. 3.2.2) . Bottom, l eft to right: (c) evolved star with 
surrounding multi-layer envelope and plume emission in the K band (Sect. |3.3| >; (d) a clumpy dust torus in 
the active galactic nuclei NGC1068 in the K band (Sect. |3.4l2) . The simulated images are presented in their 
largest field of view when available. 



The first set of images consisted of planets (top left image of Fig. [3]) creating localized 
shadows ( |Jang-Conde"lIl|2QQ9| ): (i) 10 and 20 M e planets at 8 AU, (ii) 10 and 50 M e planets 
at 4 AU, (iii) a 20 M e planet at 2 AU, and (iv) a 50 M e planet at 1 AU. All these planets 
are present in the image. The second set of images (inset of the top left image of Fig. [3]) 
were of gaps created by tidal forces between the planet and disks, similar to those shown 
in hydrodynamic simulations by |Bate e t al (2003 ). These gaps were modeled as ad-hoc 
azimuthally-symmetric perturbations. If the unperturbed surface density profile is defined 
as 2b (r) = f p(r, z) dz then the surface density profile of a disk with a gap whose center is 
located at radial distance r g of depth d and width w is 2 = 2b [1 - dexp(-(r - r g ) 2 /2w 2 )]. 
Holding 2 fixed, the vertical density and temperature structures were recalculated using the 
radiative transfer model described above. Two images of a gap at r g = 10 AU was calculated: 
(1) "big gap" scenario with d = 0.9 and w = 1.8 AU, corresponding to a 110 M e (0.35 Mj) 
planet; and (2) "small gap" scenario with d = 0.5 and w = 1.0 AU, corresponding to a 37 
M e (0.12 M 7 ) planet. 
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3.2 Granulation in late-type stars due to large convective cells 
3.2.1 Context and major issues 

Outside the spots, the solar surface is covered by bright granules with warm upflowing ma- 
terial surrounded by dark, cool inter-granular lanes where matter falls downward. Local 
multi-dimensional radiation hydrodynamics (RHD) simulations of stellar surface convec- 
tion, solving the coupled equations of hydrodynamics and non-local radiation transport in- 
cluding constant gravity, ionization effects, and realistic (binned) opacities, first successfully 
reproduced the properties of solar granulation ( Nordlund, 1982 ). This includes pattern and 
dynamics (incl. size and timescales) of the granulation, as well as center-to-limb variation, 
and line profiles ( [Dravins et al[ |1981| ). The excellent agreement between observations and 
models allowed these methods to be extended with some confidence to other stars, replac- 
ing classical ID stellar atmospheres for high-accuracy abundance determinations ( |AUende| 
|Prieto et al[[2QQT]|Caffau et al||2008) . 

The size and also the contrast of the surface structures can be directly linked to the 
efficiency of convection, to the formation of molecular lines (CO on the Sun: Wedemeyer- 



Bohm et al||2005| water on Arcturus: |Ryde et al||2002| water on Betelgeus e: |Ryde et al 



2006) , to the micro- variability (limiting the detectability of planets | Ludwig[ |2006| ) and to 



the vertical extent of "overshoot" of convective influence (Freyta g et al |1996 » into the outer 



layers (molsphere, wind: Sect. |3.3| and dust formation: Freytag and Homer 2008 ). 

Simulations of solar-type stars show the granulation properties known from the Sun, 
with a horizontal size scaled with the surface pressure scale height and slightly different 
contrast (e.g. Robinson et al, 2004). However, the situation changes for objects near the 
"granulation boundary" (A-type stars, cepheids, red supergiants): |Gray| ( |2008| ) observes 
reversed C-shaped line bisectors on Betelgeuse and relates them to granulation and giant 
convection cells accompanied by short-lived oscillations . Also Stoth ers and Leung] ( |1971| ) 
and [Schwarzschild ( 1975 ) attribute the observed photometric variability in red supergiants 
to convection cells larger than one would expect from scaled solar granulation, which is 
confirmed later by interferometric observations of global surface structures on Betelgeuse 
( |Buscher et al| [1990} [Young et al||2000| ). |Josselin and PlezlpOOT] ) find large line shifts and 
asymmetries in RSG spectra and attribute them to the existence and motion of a few large 
granules. Local RHD models of surface patches of A-type stars and cepheids (Freytag et al| 
2008 ) and global models of red supergiants (using a spherical potential and including the 
entire outer convective envelope and the atmosphere, Freytag et al 2002 ) also show similarly 
large cells and granulation with properties deviating from the solar case. 

For a large group of stars near the granulation boundary - from the main-sequence to 
red supergiants - our benchmark star, the Sun, is therefore no longer useful. There is a 
clear need for new references, and interferometric imaging will be capable of providing one 
for the important case of red supergiants, shedding light on this particular kind of surface 
granulation and its relation e.g. to the outer layers and wind formation. Recent images of 
Betelgeuse obtained with IOTA in the H band ( Haubois et al , 2009 ) display two low contrast 
spots, one of them probably marginally resolved. The spots are consistent with signatures of 
large convective cells ( Chiav assa et al[ [2009 ). 

3.2.2 Stellar surface simulation 



|Chiavassa et al (2009) computed intensity maps at different wavelengths (corresponding 
to broadband filters H and K) with the radiative transfer code OPTIM3D from snapshots 
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of the 3D hydrodynamical simulation (code CO 5 BOLD) of a red supergiant star presented 
above. They used a model with stellar parameters close to those of Betelgeuse (Levesque 
|et al[ |2005] ): a 12 M© stellar mass, a luminosity averaged over time of L = 93000 + 1300 L , 
an effective temperature of T e ff = 3490K, a radius of R = 832R , and surface gravity 
log(g) = -0.337. The numerical resolution is 235 3 with a grid spacing of 8.6R©. 

The granulation pattern has a significant impact on interferometric observables. In order 
to derive their characteristic signature, Chiavassa et al (2009) computed visibilities and clo- 
sure phases from these intensity maps. They found that convection-related surface structures 
cause fluctuations in interferometric observables which manifest themselves mostly at high 
spatial frequencies after the first null point in the visibility function. In general, the visibili- 
ties and closure phases deviate greatly from the uniform disk or limb-darkened cases, due to 
the small scale structure on the model stellar disk. These observables have been compared 
to existing Betelgeuse observations ( [Young et al] [2000| [20041 |Haubois efal|[2009| ). The 3D 
simulation gives good fits to the observed visibility points and closure phases providing a 
consistent solid detection and characterization of the granulation pattern from the optical to 
the near-infrared (Chiavassa et al, 2010| ). 

The granulation pattern can be characterized with today's interferometers by searching 
for angular or temporal visibility fluctuations or by looking at different spectral regions 
corresponding to spectral features and continuum. High-angular resolution images should 
verify or contradict the 3D simulations in terms of surface intensity contrast, granulation 
size and temporal variation. The figure at the top right side of Fig. [3] shows the image, 
computed in a narrow K band channel, that has been used in this study. 



3.3 Molecular layers around late-type stars 
3.3.1 Context and major issues 

The late stages of stellar evolution are well studied because they are the source of most of 
the dust produced in the universe ( Sedlmayr[ |1994| ). Dust is detected at large distances and 
stars are known to produce molecules which are the building blocks of dust particles. A 
molecular shell containing water vapor and carbon monoxide was first suspected around red 
supergiants by |Tsuji| |2000 ) who coined the name molsphere to describe it. This explanation 
was shown to be consistent with interferometric data collected in the near and mid-infrared 
( |Perrin et al[ |2004a[ |Weiner| |2004[ |Ohnaka| |2004| ) and subsequent works showed the list 
of molecules could be further extended to include SiO and AI2O3, the latter being a dust 
species that can survive the high temperature range of the molsphere (1500-2000 K) and 
provide nucleation sites for SiO ( |Perrin et all|2007[|Verhoelst et al||2006| ). 

Molecules in the atmosphere of Mira stars are probably lifted by large-amplitude pulsa- 
tions. This is possibly a different process from the one at play in supergiants, yet similarities 
have been found between the two classes of objects. Diameters of Mira stars are known 
to vary by factors of a few in and out of TiO bands (e.g. |Labeyrie et al[ |1977| ) as well as 
in other bands although with smaller amplitude. Measured interferometric diameters have 
been difficult to explain with models partly because of pulsations. Some dynamic models 
have been compared to interferometric data (e.g. Woodruff et al, 2004; Fedele et al[|2005] }. 
Using the Scholz and Takeda ( 1987 ) model and removing the contribution of the continuum 
in the K band, Perrin et al ( 1999) have shown a good agreement with R Leo interferometric 
data raising the hypothesis that molecules are a major contributor to the infrared surface 
brightness distribution of these objects. This was confirmed by fitting K and L band data 
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using a model incorporating a molecular shell a stellar radius above the photosphere (Men- 
ness on etal| |2002; Perr in et a!||2QQ4b|). Recent wo rks have provided more detailed studies 
of the shell (e.g. |Ohnaka et al||2QQ5| ); |Eisner et al|pQQ7] ); [Wittkowski et al|pQQ7| ); pinaka| 
et al (2009 )). In parallel models have been made more sophisticated and in particular now 
provide a physical explanation for molspheres ( Woi tke etaI| [T999). 

All studies had initially hypothesized circular symmetry and homogeneity for the layer. 
Images would yield a refinement of the structure of the layer and of the star, as well as the 
layer brightness distribution. The images of T Leporis in narrow spectral channels across the 
H and K bands taken with VLTI (Le Bou quin et al| 2009 ) have confirmed the presence of 
the circumstellar layer. The ring appears circularly uniform to within the precision of image 
reconstruction with sparse apertures, consistent with the type of shell model proposed by 
|Perrin et a l (2004b ). Images taken with IOTA of^ Cygni in the H band do not disclose such 
a bright shell, although it is detected by model-fitting to visibilities. The star itself shows 
varying asymmetries as a function of the pulsation phase ( [Lacour et al[ |2009| ). However 
insufficient dynamic range in the reconstructed image does not allow any departure from 
circular symmetry of the shell to be measured. 



3.3.2 Modeling molecular layers 

Interferometric observations (and thus also reconstructed images) of molecular layers can 
be compared either to self-consistent (hydrodynamical) atmosphere models, or to flexible 
analytical models which aim at a reproduction of the observed intensity distribution with a 
minimum of physics. 

The first approach works reasonably well for Mira stars. For example, Wittkowski et al 
(|2008|) use the mo dels by |Ireland e7all ( [2004| ), and |Woodruff et al|p009l ) use the models by 
Keller and Wood (2006 ), both on O-rich AGB stars. Theoretical preparatory work on C-rich 
AGB stars has been published by Paladin fet al| ( |2009| ), based on mode ls by | Hofner et al| 
( 2003 ), and a confrontation with observations is in progress ( |Sacutoetal[|2011| ). 

In some cases, a more ad-hoc approach is required. This usually implies a parametrized 
distribution of matter and opacity, and some radiative transfer calculation. These models 
should then be checked against thermal and chemical equilibrium. Such an approach has also 
been used with AGB star observations (e.g. |Mennesson et al[[2002||Perrin et al[[2004b|[Le 



Bouquin et al 



et al 



et al 



2009|), but it is hitherto the only option for red supergiant stars (e.g. Verhoelst 



2006; Ohnaka et al, 2008 ), or very asymmetrical circumstellar environments (Kervella 



2009) 

The synthetic images used in this work consist of an RSG photosphere model, computed 
with the marcs code ( [Gustafsson et al[[T975[ and further updates), surrounded by two layers 
of molecular material: a layer of CO is located at 1 . 15 R*, with a temperature of 2350 K, and 
a column density of 10 23 cm~ 2 . The second layer contains 10 21 cm~ 2 of H 2 and is located 
at 1.7 R*, with a temperature of 1850K. Around this object, we placed a dust shell which 
reproduces the observed silicate features in QfOrionis (M ~ 3 x 10~ 8 M© yr _1 ). This shell 
starts at 13 R*, and has a typical r~ 2 outflow density profile. To mimic the recent results of 
|Kervella et al| ( [2009] ), we included a plume as a characteristic of inhomogeneous mass loss, 
possibly related to convection, as recently observed in Betelgeuse. The chosen plume shape 
is a simplification of a structure originating at the photosphere, experiencing dilution as it 
travels outwards. It represents a local increase of 50% of the circumstellar intensity. Images 
were computed across the ^Tband, but only the one near 2.4 yum is presented here (bottom 
left side of Fig. [3] and zoomed image at top centre of Fig.|6|, as it shows both the CO and 
the H 2 layer. 
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3.4 Dusty tori in active galactic nuclei 
3.4.1 Context and major issues 



Although there are numerous classes of active galactic nuclei (AGN), they have been uni- 
fied into a single scheme (e.g., Antonucci, 1993 , 2002; Urry and Padovani, 1995 ). The basic 
premise of unification is that every AGN is intrinsically the same object: a super-massive 
black hole whose activity is powered by accretion through a surrounding disk. This central 
engine is further surrounded by a dusty toroidal structure, which provides anisotropic ob- 
scuration of the central region. The observed diversity is simply the result of viewing this 
axisymmetric geometry from different angles. Directions with clear sight of the central en- 
gine yield "type 1" objects, where clouds in close proximity to the black hole, moving at 
high velocities because of the strong central gravity and ionized by the strong central contin- 
uum radiation, give rise to broad line emission in high ionization lines. AGN whose central 
engine and broad-line region (BLR) are blocked from view are "type 2". 

The primary evidence for the torus arose from spectropolarimetric observations of type 
2 sources, which reveal hidden type-1 emission via reflection off material situated above the 
torus opening. While compelling, this evidence is only indirect in that it involves obscu- 
ration, not direct emission by the torus itself. Direct evidence comes from IR observations 
— an obscuring dusty torus should reradiate in the IR the fraction of nuclear luminosity 
it absorbs, and the continua from most AGN indeed show significant IR emission. Recent 
progress in IR high-angular resolution techniques achieved direct imaging (e.g.|Weigelt et al| 
2004 ) and interferometric measurements ( [Jaffe et al]|2004|[Mei senheimer et al , 2 007 [ Tn£\ 



tram et al, 2007 ) of the torus in a few AGN. These observations show that the torus is rather 



compact. The size scale is set by the dust sublimation radius R& 0.4 (L/10 45 erg s _1 ) 1/2 pc, 
which determines the torus inner radius. All observations are consistent with R 0VLt /R^ no 
larger than ~ 20-30, and perhaps even only ~ 5-10. The observations also indicate close 
proximity between dusty regions with widely different temperatures. These findings can 
only be understood if the dust distribution is clumpy because in that case the dust temper- 
ature is much higher on the illuminated face of an optically thick cloud than any other part 
of its surface, resulting in co-location of different dust temperatures. The torus dumpiness, 
predicted early on by Krolik and Begelman ( 1988 ), naturally explains many features of the 
IR observations that cannot be understood in smooth-density models (Nenk ova et al| [2002; 
|Elitzur[ [20081 |Nenkova et al||2008a|b) . 



3.4.2 Clumpy tori simulations 

We simulated model images based on 3D clumpy torus models ( Honig et al , 2006 ; Honig 
|and K ishimoto , 2010). The main parameters of the models are (1) the radial dust distribution, 
parametrized as a power-law n(r) oc r°, (2) the mean number of clouds No along an equatorial 
line-of-sight, and (3) the half-opening angle #o (defined as cr in a Gaussian distribution 
along the altitudinal direction). The inner radius of the torus R& is set by dust sublimation 
at temperature T su b = 1500 K and scales with the luminosity of the AGN, R& oc L 1/2 . While 
the model spectral energy distribution (SED) has a negligible dependence on the actual size 
of the dust clouds (see Honig and Kishi moto[|2010[ Sect. 2.5), the clouds might show up 
in the interferometric signal as small-scale emission or obscuration regions. The cloud sizes 
are parametrized by a radial power-law. 

Starting from pre-calculated databases, clouds are randomly distributed within the torus. 
Each cloud is associated with a corresponding cloud in the database of directly- or indirectly- 
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Table 2 Imaging interferometers in use or in construction. 



Facility 


Wavelength 


Numbers of 


Aperture 


Baseline 




(microns) 


apertures 


(m) 


(m) 


NOI 


0.45 - 8.5 


6 


0.12 


2 - 99(437) 


CHARA 


0.5 - 2.5 


6 


1 


50 - 330 


VLTI 


1 - 13 


4 


8.2 


40-130 




1 - 13 


4 


1.8 


8-200 


MROP 


1-2.5 


10 


1.4 


7.5 - 350 



* in construction 



illuminated clouds, the cloud image is placed in the model space, and images at all spectral 
channels are computed using the obscuration along the observer line-of-sight. The simula- 
tions are performed for a torus which is seen edge-on by the observer, as we would expect 
in a type 2 AGN, located at the distance of NGC 1068, the only AGN with published inter- 
ferometric data in both near- and mid-IR ( jWeigelt et al) [20041 [wTttkowski et a!) [20041 |Jaffe| 
|et al , 2004; R aban et al||2009| ) which is close enough so that dumpiness in the torus can be 
observed (see the right image of Fig.[3]to get a wide view of the simulated AGN and the top 
right image of Fig. [6] to get a zoom on the central part). For that simulation we did not take 
into account the incoherent flux contribution from the surrounding galaxy which will bias 
the measurements. We discuss this issue further in Section [5X4] 



4 Available imaging facilities 

NOI, CHARA, VLTI, and MROI are the present and near-future interferometers which can 
perform image synthesis (see Table [2}. We used their configurations to produce realistic 
simulated reconstructed images for the astrophysical topics described in Sect. [3] We describe 
below the main characteristics of these instruments and their most recent imaging results. 



4.1 CHARA 

The Center for High- Angular Resolution Astronomy (CHARA) Array is a 6 telescope inter- 
ferometer operated by Georgia State University on the Mt. Wilson site outside Los Angeles, 
California (ten Brummelaar et al, 2005). CHARA has the longest baselines of any current 
optical interferometer (maximum 330m baseline), allowing sub-milli-arcsecond angular res- 
olution in the visible and infrared. While most of the 30 CHARA refereed papers to date 
present infrared observations using only single baseline data (e.g., McAli ster et al[ [2005) , 
we focus here on progress in imaging which requires at least 3 telescopes to be used. 

Currently, the only combiner at CHARA to have published imaging results is the Michi- 
gan Infrared Combiner (MIRC, |Monnier et alj|2006| ). MIRC now combines up to 6 CHARA 
telescopes simultaneously and can produce snapshot images of stellar surfaces and interact- 
ing binaries, as long as the surface features are not too complicated. Mo nnier et al| ( |2007| ) 
presented a detailed image of the star Altair, the first resolved image of a main-sequence 
star, confirming the strong gravity darkening and oblateness expected for this rapid rota- 
tor. The rapid rotators Alderamin and Rasalhague have now also been imaged and modeled 
by CHARA ( [Zhao et al| |2009). Another breakthrough made possible by imaging was the 
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first "movie" of the interacting binary ft Lyrae ( [Zhao et ai| |2008), imaging the pair of stars 
through half of the ~13 day orbital period. More recently impressive images of the well 
known eclipsing system e Aurigae revealed the transiting shadow of a flat disk ( Klo ppen-| 
borg et al , 2010 ). New imaging targets include the gas disks of Be stars and rotating spots 
on active stars. 

Imaging efforts at CHARA are rapidly intensifying with the recent commissioning of 
the PAVO and VEGA combiners, both visible combiners capable of closure phase mea- 
surements and interferometric imaging. A 6-telescope fringe tracker (CHAMP) has been 
delivered in Summer 2009 and the 2-telescope CLASSIC combiner has been upgraded to 
use three telescopes (CLIMB). MIRC is starting to use all six CHARA telescopes simulta- 
neously allowing the imaging of planet-forming disks in young stellar objects for the first 
time. 



4.2 VLTI 

The Very Large Telescope Interferometer (VLTI,[Haguenauer et al, 2008) is located at Cerro 
Paranal, Chile, and is operated by the European Southern Observatory. The infrastructure 
is composed of four relocatable 1.8m Auxiliary Telescopes (ATs) equipped with visible- 
wavelength tip/tilt correction, four fixed 8.2m Unit Telescopes (UTs) equipped with visible- 
wavelength adaptive optics, and 6 delay lines of 100 m stroke. The minimum baseline is 8m 
while the maximum one is about 200m. The focal laboratory is equipped with a JHK-band 
3-beam combiner AMBER ( [Petrov et~al| |2007| ) allowing precise measurements of closure- 
phases, and with an Af-band 2-beam combiner MIDI (Leinert et al, 2003 ). The VLTI instru- 
ments are opened to the astronomical community thanks to dedicated observatory support, 
leading to a total of more than 191 refereed papers in 8 years (2004-201 1). 



The first reconstructed image produced by the VLTI is of the binary 6 l Ori C ( Krausetal 
2009), using three configurations of 3 ATs. The spectral coverage of AMBER across the K- 
band was used to enhance the coverage of the (u, v) plane. However because observations 
were not close enough in time, the authors had to artificially "rotate" the (u, v) plane between 
each configuration to account for the binary rotation. Almost at the same time, Le Bouquin 
|et al| ( |2009] } reconstructed a spectro-image of the Mira variable T Lep at a spectral resolution 
R ~ 35 using AMBER and four 3 -AT configurations. Using the MIDI instrument with the 6 
possible UT baselines, Raban et al (2009) reconstructed an image of the torus at the heart of 
NGC 1068. The image is point- symmetric because of the lack of phase information. Other 
imaging works are in progress concerning B[e] stars (Mi llour et al[ |2009| ), young stellar 
objects and supergiants. 

ESO is currently commissioning PRIMA (Delplancke et al, 2006), an instrument which 
will be able to simultaneously observe two stars separated by up to an arc minute, using 
2 telescopes. Beyond the astrometric capability which is the primary focus of PRIMA, it 
will allow differential phase measurement and possibly an improvement in absolute phase 
measurements for imaging. The fringe sensors that come with it are also expected to be 
more sensitive. Since October 2011 a visitor instrument PIONIER allows four VLTI tele- 
scopes (ATs or UTs) to be combined and has provided first images of interacting binaries 
( [Le Bouquin eTal) [2QTT] |Blind et al[ |2011b| |Berger et al| |20T0] ). Around 2013, the VLTI 
is expected to receive the first second-generation general user instruments making use of 4 



telescopes simultaneously. These are Matisse working in the L, M, N bands (Lopez et al 



2008 
2008 



and Gravity with dual beam operation and an astrometric capability (Eisenhauer et al 
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Later, it is planned that the VSI instrument will eventually combine 6 or even 8 tele- 
scopes simultaneously, providing realistic snapshot imaging capability ( Malbet ^tal] |2008| ) . 
In parallel, phase A studies of a second generation fringe-tracking instrument, capable of 
cophasing an array of 4 and potential more telescopes are being conducted. This instrument 
will be essential to improve the sensitivity of second generation science combiners ( Meisner 
[eTa!l|2QTQl [Blind et al||20lTi| ). 

4.3 MROI 

Located in New Mexico (USA), the Magdalena Ridge Observatory Interferometer (MROI, 
Creech-Eakm an et"al] |2008| is a visible and infrared telescope array under construction. 
MROI will comprise ten 1.4-m diameter unit telescopes laid out in a Y-shaped array. The 
telescopes will be relocated approximately three times per year. Four scaled array configu- 
rations will be available, the most compact configuration giving shortest baselines of 7.8 m 
for comparison with single-dish measurements, and the most extended configuration having 
a longest baseline of 340 m. 

MROI has been designed for model-independent imaging of faint targets with a layout 
optimized for baseline bootstrapping employing equal spacings between adjacent telescopes 
on each array arm, essentially 100% sky coverage down to 60 deg from zenith, a beam train 
that minimizes the number of optical surfaces encountered by the light beams, and a group- 
delay fringe tracking capability, with limiting point-source sensitivity of H - 14. 

First light is expected in 2012, with first fringes in 2013. Phase A deployment (6 tele- 
scopes and the JHK science capability) could be completed as soon as 2015. Phase B (10 
telescopes and optical science) will begin when more funding for unit telescopes, associated 
delay lines and an optical beam combiner is obtained. 

4.4 NOI 

NOI (ex NPOQ) has pioneered aperture synthesis in the optical domain and is still oper- 
ational today with active plans for future developments ( Armstr ong et al[ 1 1998b} |Hutter| 
|et al| |2008| ). The array has two subarrays dedicated to imaging and astrometry. While the 
imaging array uses 12 cm beam-stopped relocatable siderostats it will eventually use six 
siderostats of 35 cm aperture. The astrometric subarray is composed of six 50 cm siderostats 
with shorter baselines (19 to 38 m). The focal instrumentation operates in the 450 nm to 
850 nm band. The current imaging capability is obtained through the combination of two 
"imaging" siderostats with four of the astrometric ones. NPOI offers the unique feature 
of being able to maintain active group delay fringe tracking for its 6 siderostat operation 
therefore offering a unique snapshot imaging capability. The current maximum operational 
baseline for imaging is 99 m but the array will have the capability to span baselines from 
2 to 437 meters. The integration of the former Keck-interferometer 1.8 m telescope out- 
riggers to the imaging array is under study. In parallel a study to install 1.4 m composite 
telescopes equipped with adaptive optics is underway with the first step being construction 
of a prototype. Since 1998 NPOI has been particularly successful at improving its automa- 
tion (in particular in the field of array -phasing), studying binary systems of various types, 
and demonstrating a spectral phase referencing capability (e.g. |Benson et al[ |1997| [Zavala 
[e7a!l[20T0irSchmitt et al||2009] ). 



6 Navy Prototype Optical Interferometer 
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Fig. 4 From an image to visibilities: Modeling interferometric observables with ASPRO. 
5 Results from image reconstruction simulations 

In this section we first describe the method for generating realistic interferometric data from 
the MROI and VLTI, utilizing the telescope relocation capability of these facilities. Then we 
show the images that were reconstructed by our colleagues without prior knowledge of the 
scientific content. Finally we compare the reconstructed images with the original ones and 
discuss our findings. 



5.1 Simulating realistic interferometry data 

The simulations used here have been made with the ASPRO package ( [Duvert et al[ |2002| ). 
ASPRO is a general-purpose observation preparation tool developed by the Jean-Marie Mar- 
iotti Cente^jfor optical interferometers, with a particular emphasis on the VLTI instruments 
(Sect. |4.2| . ASPRO and its little brother ASPRO 2 are available as a web-based utility, with 
a client-server architecture and a Java user interface, at http://jmmc.fr, and also as a 
standalone program. ASPRO simulates the observation by an interferometer of an astronomi- 
cal object, at one or several times, and delivers simulated interferometric observables in the 
OIFITS format f Pauls et al[|2005] }. 

The interferometric chain is modeled as the combination of an interferometer infrastruc- 
ture and focal instruments. The interferometer infrastructure comprises the telescopes, delay 
lines, tip-tilt correctors, adaptive optics, and fringe trackers. It adds geometrical require- 
ments such as the positions and sizes of the telescopes apertures at the time of observation 
with respect to the position of the science object and the geometrical delays thus induced 
between each pair of telescopes. In addition, it includes environmental constraints such as 
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the atmospheric seeing, the different horizons seen by telescopes and technical limitations 
(limits on the delay-line strokes, flux dependence on the performance of active optical el- 
ements...). For each interferometer simulated in ASPRO these requirements are known and 
tabulated; only the atmospheric seeing (and its associated coherence time) is a variable factor 
that can be adjusted by the user. 

The focal instrument is described by its spectral capability (bandwidth, resolution, num- 
ber of spectral channels, central wavelength, photometric band...) and its sensitivity and 
throughput. The former are, for most instruments, fixed as a list of instrumental "modes". 
For example, for AMBER, there are 18 such modes, one giving simultaneous /, H, K obser- 
vations at 35 resolution, 3 for medium (1500) resolution in H or K, and 14 for high (12000) 
resolution in K. The sensitivity and throughput, i.e., the limiting magnitudes, the visibility 
losses, the biases and the variance of the values of the various interferometric observables 
are computed with a dedicated noise model (Appendix |A). 

The science object is described by its position on the sky and its angular and spectral en- 
ergy distributions. Those can be either a series of images, one for each wavelength, grouped 
in a data cube, or a combination of simple analytic models (binary, uniform disk, etc.). 

5.2 Simulation methodology 

The simulation of observations, depicted in Fig. [4] consists of: 

1 . building the list of (u, v) positions of the baselines at each wavelength of the required 
focal instrument mode and at the various times (for example, an observation every 20 
minutes for 6 hours around the meridian) intended by the observer; 

2. computing the Fourier Transform (FT) of the flux spatial distribution of the object at 
each wavelength at each (u, v) point defined above (in the present case, where the user 
provides a model of the object as a cube of images); 

3. from the complex visibilities obtained, computing all the interferometric observables, 
V 2 , amplitude and phase of the bispectrum, and, if the instrument has spectral capabili- 
ties, the amplitude and phase of the differential visibility; 

4. for each of these interferometric observables, applying the noise model (Appendix |A|) 
and a lower- limit accuracy threshold to get the expected variance of each value; 

5. randomize the interferometric observables according to their expected variance; 

6. export these simulated observations in the OIFITS format. 

The (u, v) coverage for both MROI and VLTI was derived from the allowed 6 telescope 
configurations (Tableland Fig. [5]). In the case of MROI this was obtained by switching 
rapidly 6 beams in front of a four beam combiner. This explains why threexfour telescope 
combinations are used to get all 15 baselines. 

Undoubtedly, because we do not take into account a detailed model of fringe tracking 
performances (outside the scope of such a paper), we set ourselves on the optimistic side 
in terms of calibration biases. We can therefore consider that the images obtained here are 
the best one can expect from the interferometers and that actual data may contain additional 
artifacts. 

By simulating the expected results of the actual observation and providing OIFITS files, 
ASPRO allows assessment of not only the feasibility of the observations but also the sci- 
ence potential of the planned observations. The OIFITS file obtained can be used to check 
whether image reconstruction techniques retrieve the desired information and to compare 
the results to what a model-fitting software can achieve. The interferometric observables 



Imaging the heart of astrophysical objects with optical long-baseline interferometry 



23 



Array 


Configuration 


MROI 


N2 W2 W4 S2 




N4 W2 W4 S4 




N2 N4 S2 S4 




N4 W4 W6 S4 




N6 W4 W6 S6 




N4 N6 S4 S6 


VLTI 


A0B1 D2G1 J2 MO 




Al B3 D2 Gl HO 11 




UT1-UT2-UT3-UT4 



Table 3 Array configurations for both MROI and VLTI used in the simulations 
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Fig. 5 The (u, v) plane coverage used for the image reconstructions. Top from left to right: Matisse/VLTI-UT 
coverage for the protoplanetary disk model; VSI/VLTI-AT coverage for the late-type supergiant star model for 
Betelgeuse; VSI/VLTI-AT coverage for the molecular layer model around a late type star; I4/MROI coverage 
for the dusty tori model for the active galactic nuclei NGC 1068. 

Table 4 Description of the parameters of the different science cases 
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Table 5 Number of squared visibilities (V 2 ), closure phases (CP) and their corresponding average errors over 
the N V 2,N CP 2 measurements. To simulate a potential calibration limit in the minimum visibility and closure 
phase that could be measured we introduced a threshold of 0.002 in visibility and 0.1 degree in closure phase. 
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can also be compared to the model observables in terms of chi-square distance as a first 
estimate on the feasibility of the planned observation. Table [5] contains a summary of the 
computed errors on the visibility and closure phase for the different input images. 

5.3 Results 

We have performed a series of simulations following the astrophysical prescriptions given 
in Sect. [3] and using the different instrument configurations described above. We focus our 
attention on the four cases listed in Table]?] since we chose not to explore all the possibilities, 
which would have led to more than several tens of results. We have selected results obtained 
with the instruments VSI/VLTI and I4/MROI in the K band and Matisse/VLTI in the N band. 
The (w, v) coverage for each configuration is displayed in Fig. [5] The configurations were 
chosen visually to fill the (u, v) plane as uniformly as possible and therefore ensure the best 
spatial frequency sampling. For the MATISSE case a fixed UT configuration was used. For 
the VSI reconstruction two 6 AT telescope configurations were chosen with the supergiant 
one pushed to the maximum possible hour angle range (e.g -5h,5h). The MROI coverage 
corresponded to one of the rare AGN with sufficient angular extension to permit imaging 
(NGC 1068). 

The Planet signatures in protoplanetary disks case (hereafter called the protoplanetary 
disk case) was the only one for which we were not able to retrieve a viable reconstructed im- 
age. The conclusion applies both to the embedded-planets image and the disk+gap images. 
The explanation for this lack of success is that the high contrast between the planet induced 
gap or shadows and the disk/photosphere emission clearly exceeded the dynamical range 
achievable with such a technique. None of the cases computed in Sect. |3. 1.2] was successful 
in imaging features beyond the central star, using BSMEM or MACIM. However, the ex- 
ample chosen here is voluntarily an extreme case of imaging. Long -baseline interferometry 
has already begun to explore imaging of protoplanetary environments (e.g. |Benisty et al| 
|2011| ) and it is expected that MATISSE will have the ideal combination of sensitivity and 
angular resolution to reveal the inner morphology of dust emitting inner disks. Possible im- 
provements to the method would be to introduce a multi- wavelength approach. The location 
of the gap being independent of wavelength continuum observations at various bands might 
provide additional hints to constrain the gap positioning. We note that recent gaps in tran- 
sition disks were indirectly revealed by visibility modeling (Benisty et al, 2010; Olofsson 
|eTall|2QTTT >. 

For the remaining cases the results are summarized in Fig. 6 and Fig. 7 with the simu- 
lated images (top row), the simulated images convolved by a Gaussian beam at the maximum 
interferometer angular resolution (central row), and finally the reconstructed images (bottom 
row). Figures [8] and [9] show a contour comparison between model and reconstructed image. 
These figures allow us to carry out the first qualitative comparison between the reconstructed 
and original images at the proper resolution. For the Molecular layer around late-type stars 
case (hereafter called late-type envelope case) and the Dusty torus around AGN case (here- 
after called the AGN case), note that the images have been magnified in order to display the 
central part of the images. In these cases, the image reconstruction was not able to detect the 
less luminous structures present in a wider field (Fig. [3]). 

A first visual inspection shows that, 

1 . for all cases, all the features with flux greater or equal to approximately 1/20 of the peak 
intensity in the image are roughly correctly positioned and have the correct shape and 
size; 
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Fig. 6 Resulting images from the image reconstruction (bottom row) of the supergiant (left) and evolved 
star (right) cases compared to the initial images resulting from our modeling (top row) and compared also 
to the image that would be obtained by a single telescope whose diameter is the maximum interferometer 
resolution (central row). The astrophysical cases are on the left column a supergiant star and on the right 
column a molecular layer around late type star. The color codes correspond to logarithmic scales and the 
image have been zoomed-in with respect to figure[3]to concentrate on the meaningful reconstructed details. 
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Fig. 7 Resulting images from the image reconstruction (bottom row) of the AGN images compared to the ini- 
tial images resulting from our modeling (top row) and compared also to the image that would be obtained by 
a single telescope whose diameter is the maximum interferometer resolution (central row). The astrophysical 
object is here a dusty tori in active galactic nuclei. The color codes correspond to logarithmic scales. 



Imaging the heart of astrophysical objects with optical long-baseline interferometry 



27 



2. there is evidence that the image reconstruction process brings slightly better resolution 
than the fringe spacing A/B which is due to the super resolution allowed by the regular- 
ization process and the interferometer response function (this result is best illustrated in 
the AGN case in Fig. [7]); 

3. in the supergiant case, the most prominent "granular" atmospheric features are retrieved 
but the dark thin intergranular lanes of the underlying large-cell convective structure, 
visible in the model image, are not; 

4. in the late-type envelope case the basis of the plume on top center of the image is 
marginally detected; 

5. in the AGN case, the complexity of the scene is obvious and the emission morphology 
is globally well-reproduced in the reconstructed image. However as we will see later the 
quantitative interpretation of such an image requires caution. 

A proper quantitative comparison cannot rely solely on this visual inspection. For the 
discussion below, we used one possible image reconstruction quality estimator, called fi- 
delity, described in Appendix [B] and inspired from ALMA work. There are also valid met- 
rics like the variance of the difference between the two images, reconstructed and original 
smoothed at the interferometer resolution (Renard et al , 2011) or the figure of merit pro- 
posed by |Cotton et~al| ( |2008| ) but there is no single metric with proven superiority over all 
other methods. 



6 Discussion 

Except for the protoplanetary disk case, optical interferometric imaging is quite successful in 
retrieving images similar to the true images. In this section, we discuss how new and useful 
astrophysical information can be extracted from interferometric imaging, the similarities and 
differences with radio-interferometry synthesis imaging, and how to further improve image 
reconstruction. We conclude with a discussion of how to improve the current facilities to 
obtain the best information possible using aperture synthesis imaging. 

6.1 Getting new and useful astrophysical information 
6.1.1 Planetary signatures in protoplanetary disks 

High-resolution imaging of protoplanetary disks has the potential to reveal where and when 
planets form as can be seen on the left image of Fig. [3] Growing cores of giant planets 
gravitationally and tidally perturb the protoplanetary disks in which they form to create 
indicators of their presence on scales much larger than the planets themselves, including both 
localized dimples and azimuthally symmetric gaps. Shadows and brightenings in scattered 
light images from these perturbations will reveal their presence. 

The difficulty in detecting these perturbations using optical interferometry is that the 
stellar and disk brightness overwhelms the emission from the features in the disk. At mil- 
limeter wavelengths, the stellar emission is minimal, so the dynamic range is less of an issue 
although these wavelengths probe a different regime of the disk as well: the disk becomes 
optically thin and the images represent vertically integrated thermal emission from the disk. 
Because of the nature of shadowing and illumination at the surface of the disk, the contrast 
of the shadows and brightenings is higher at visible and infrared wavelengths as compared 
to radio wavelengths. 
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Fig. 8 Contour plots of both original (top row) and reconstructed image (bottom row) for the supergiant case 
(left) and evolved star with plume (right). Contour levels are expressed in fraction of maximum intensity. 
Contours of bottom row correspond to the model image for better comparison.Levels: 0.9,0.8,0.7,0.5,0.3,0. 1); 
evolved star with plume (levels: 0.9,0.5,0.3,0.21,0.1). 
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Fig. 9 Contour plots of both original (top row) and reconstructed image (bottom row) in the AGN science 
case. Contour levels are expressed in fraction of maximum intensity. Contours of bottom row correspond to 
the model image for better comparison.Levels: 0.9,0.5,0.1,0.01,0.001). 
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However, these objects can still bring some interesting surprises. In the past, it was 
believed that disks around young stars would be only marginally resolved by long-baseline 
optical interferometers in the infrared, but in fact the structures which have been measured 
are indeed located at larger scales than expected (see review bypillan-Ga bet et al| [2007] ). 
There might be several other signatures of planets and their formation in the potential planet 
forming region of circumstellar disks (inner 0.1" in nearby star-forming regions): overdense 
vortices ( [Wolf and Klahr] [2QQ2|); gaps ([Wolf et al| |2QQ2} jVarniere et al| [2004]) ; accreting 
giant planets ( |Wolf and D V Angelo||2005[|Hofmann et al| |2006[|Woff et al(|2007| > and other 
prominent features. Indeed Malbet etld^ 2005| )"and Millan-Gabet et al (2006 ) have reported 
oscillations in the visibilities that could trace the presence of ofT-centered hot spots in the 
disk without identifying its nature. If confirmed, imaging of these relatively bright spots 
should be relatively easy since the dynamical range quoted in these works ranges between 15 
and 50. The failure of our tests reveals that the image reconst ructions are limite d to dynamic 
ranges significantly smaller than 10 3 . The reader is referred to Wolf et al ( 2012 ) for a detailed 
review of the potential of long-baseline interferometers for the study of circumstellar disks 
and planets. 



6.1.2 Granulation in late-type stars due to large convective cells 

The image we present here was computed in the K band. Chiavassa et al (2009 ) show that 
in this spectral region the continuum-forming layers are more visible and the granulation 
pattern is characterized by large-scale granules of about 400-500 R Q (^60% of the stellar 
radius in Fig. 4 of Chiavassa et al, 2009 ). On top of these cells, there is small scale structure 
with characteristic size of 50-100 R , ^10-15% of the stellar radius. This pattern is clearly 
visible also in the top left panel of Fig. [6] 

There are three main characteristics of the granulation pattern that we may try to recover 
from image reconstruction: 

- the size of the granulation cells from the large scale granule to the smaller scale struc- 
tures; 

- the intensity contrast; 

- the timescale of the convective motions. 

Images of stellar surfaces at various wavelengths will allow these three questions to be 
answered, and will therefore be a very valuable tool to (in) validate our simulations. This is 
crucial in order to understand the dynamics of the atmosphere of red supergiant stars, and 
its importance in, e.g., the still unknown process at the origin of mass-loss ( Josseli n~and| 
|Plez[|2007| ), or in order to determine accurate abundances from spectra using proper model 
atmospheres. 

We compared the original image, smoothed to the interferometer resolution (top left 
panel of Fig. [6), to the reconstructed image (bottom left panel). The reconstructed image 
shows medium sized granules, similar to the original smoothed image. There is, however, 
much more information in the original image, and key features like the large central granule 
are very hard, if not impossible to see in the smoothed and reconstructed image. This is due 
to the fact that the emission coming from the inter-granular lanes is too weak compared to 
the surrounding regions. 

We also compared the intensity contrast of the smoothed and reconstructed images by 
using normalized intensities. For this purpose we normalized the intensity to the peak inten- 
sity. The top images of Figures [5] and [9] show for each science case a contour display of the 
initial images, the bottom ones correspond to contour displays of the reconstructed images. 
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Comparison can be quantified by measuring the fidelity of the image (Appendix |5]). We can 
see that the fidelity of the reconstructed pixels with respect to the original ones can be fairly 
high (> 10) on a significant fraction of pixels (TableJS) but this is unrelated to the interesting 
stellar surface structures. 

Several issues need to be discussed: 

- In the K band, it is possible to recover small-to-medium scale granules, while the infor- 
mation about large convective cells is lost. While near-infrared wavelengths have been 
used in this study, one could consider the interest of studying Red Supergiants in the 
visible. In the visible, the angular resolution is higher (e.g almost 4 times between the V 
band than in the K band) and the intensity contrast of granulation is higher. At shorter 
wavelengths, corresponding to strong molecular absorption, e.g. TiO bands, synthetic 
intensity maps show more details and much larger contrast between bright and dark ar- 
eas than in the K band, and larger inter-granular lanes. While pioneering imaging results 
were carried out using visibile instruments (Baldwin et al 1996; Benson et al 1997; 
| Young et~aH|2000) near infrared wavelengths have been the dominant source of obser- 
vations because of their lowered sensitivity to atmospheric conditions. However recent 
developments at CHARA with the VEGA instrument have re-opened the way to such 
projects (e.g. jMourard et al[[2QQ9| ). 

- Test image reconstructions should be made with interferometers containing more tele- 
scopes, or larger baselines, in order to investigate the best affordable way to recover the 
granulation properties outlined above. As is suggested from comparing simulated im- 
ages in the visible and near-IR, the right choice of the wavelength may well be a key, 
despite technical difficulties growing at shorter wavelength. 

- Simultaneously to the measurement of the granulation size at different wavelengths, the 
temporal evolution of the convective pattern should also be compared to the simulations. 
Extracting typical convection turnover times from series of images should be possible 
for the medium scales, as soon as the image reconstruction can be trusted, which seems 
to be the case, even with a small number of telescopes (Fig.|6|. 

6.1.3 Molecular layers around late-type stars 

On a qualitative level, the image reconstruction was successful at detecting the outer molec- 
ular shell as a ring-like structure. The base of the plume is also visible, and the wavelength- 
dependent size of the central object corresponds to the effect of the inner molecular layer. 
Pixel fidelity computed in Appendix [B] shows indeed that the stellar- surface pixels are re- 
trieved with a good fidelity (on average 10, while the surrounding layer is less well deter- 
mined (average fidelity « 5). Our simulation did not include any dumpiness in the circum- 
stellar environment, but the present results indicate that isolated and/or large enough clumps 
with a contrast of / c ium P /^background > 0.1 would be detectable. However, the fully resolved 
extended dust structure could not be detected with the reconstruction algorithms, and the 
reconstructed limb darkening of the central object is much stronger than that in the input 
image, making diameter estimates for the central star unreliable and necessitating rather de- 
tailed visibility model fitting. Although this has not been investigated at length, we think this 
is due to the lack of short baselines. Using an image of the star obtained at the diffraction 
limit of a large telescope when performing the image reconstruction would probably solve 
this issue like single dish observation in radio astronomy. The latter problem is partly due 
to an insufficient angular resolution in this particular simulation: |Lacour et al| ( [2008] ) show 
that the limb darkening can be reconstructed if the data sample the 2nd and 3rd lobe of the 
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visibility curve, and a limb-darkened prior is chosen for the image reconstruction. Given the 
highly variable nature of many evolved stars, it is important that real observations come as 
close as possible to snapshot images. 

An evolved star and its circumstellar environments span a huge range in sizes, from 
convective cells on the surface to the dusty shell at hundreds of stellar radii. Although these 
different components lead to comparable total flux contributions at IR wavelengths, the dif- 
ference in size of the radiating surfaces immediately indicates a huge range in intensities. In 
fact, in our simulated image, the contrast between central star and inner dust shell intensity 
is /★//dustsheii ~ 10 5 , which is clearly beyond what can be retrieved from the reconstructed 
image. 



6.1.4 Dusty torus in an active galactic nucleus 

One of the most interesting tasks for interferometric imaging of AGN dust tori is the char- 
acterization of the substructure. At the moment, the actual size of the dust clouds and their 
properties is closely constrained by SED fitting (Sect. |3.4| . However, Wittkowski et al| ( [2QQ4] ) 
reported on ^-band long-baseline VLTI/VINCI interferometry of NGC 1068 where they 
found that ~40% of the emission at 46 m baseline is not coming from the larger-scale emis- 
sion region but confined to scales <5 mas. Without knowledge of the real distribution of the 
small-scale flux, this can be interpreted either as emission from the dust- sublimation region 
or accretion disk, or emission from dust clumps (or "clumping of clumps") seen through the 
obscured screen of other clouds (Honig et al , 2008 ). 

The reconstructed model images of the ^-band emission of NGC 1068 show several 
small emission regions in the central region of the torus. They reflect the parts of the emis- 
sion from the inner torus region seen through holes in the distribution of clouds at larger 
radii. The total emission contained in these regions and their sizes are consistent with the 
constraints from VINCI on NGC 1068. In this way, interferometric imaging allows us to 
directly detect substructure in the torus and even give an estimate of its contribution to the 
overall flux of the torus. Interestingly, the analysis of our model images shows that the holes 
in the obscuring screen of more distant clouds have similar sizes as the typical cloud sizes. 
If that holds for other model parameters as well, measuring the sizes of the small-scale 
emission regions might even provide constraints for the typical size of dust clouds in the 
torus. 

When trying to evaluate the fidelity of the pixel intensities with respect to the model 
images (Appendix |5]), the fidelity is generally low and while the global morphology of the 
emitting regions is nicely restored, the intensity of the emitting regions is largely deviating 
from the original one. This was expected as this image was by far the most complex of our 
set. Improving intensity accuracy is essential in order to exploit color maps (e.g H-K, K-L 
or N - H) to bring constraints on dust obscuration and the temperature distribution. 

An additional issue that was not taken into account is the probable contribution of in- 
coherent flux by the host galaxy. As recognized by Kishim oTo et al|p0TT] ), particular care 
will have to be taken in calibrating this flux through SED fitting and, when possible, through 
single dish measurements such as those provided by aperture masking (e.g. Mon nier et al| 

|2004] }. 
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6.2 Similarities and differences with mm- wavelength synthesis imaging 

While early VLBI, with its poorly calibrated phases, confronted difficulties very similar 
to optical aperture synthesis, we thought it would be interesting to compare optical long 
baseline interferometry with the millimeter one. mm-interferometry development offers an 
interesting operational parallel to what could be the future of optical interferometry. 

Interferometry and aperture synthesis in the radio domain has been routinely achieved 
since the 1970s. The last 20 years have witnessed maturation in the millimeter domain 
as well, from the construction of large facilities such as the IRAM Plateau de Bure in- 
terferometer (six 15-m antennas, Guillo teau et al[ p~992| ) to the Atacama Large Millime- 
ter/submillimeter Array (ALMA, sixty-four 12-m antennas). Compared to the infrared do- 
main, two major differences make millimeter interferometry easier to implement: (1) the 
coherence time of the atmosphere is longer, making it possible to reference the phase of the 
visibility to avoid relying on phase closure only, and (2) the millimeter receivers are directly 
sensitive to the incoming electric field, which is then down-converted to a frequency that can 
be processed with electronic equipments; as a consequence, the delay between two antennas 
can be corrected for after detection of each signal (before any interference) and the signal can 
be amplified to allow to combine many baselines simultaneously. Unlike in optical interfer- 
ometry the phase calibration can be done by observing, non simultaneously, a point-source 
calibrator (within a much larger isoplanatic patch) within one coherence time. Therefore 
the complex visibilities, including both amplitudes and phases, can be directly computed. 
Images are then obtained by Fourier transforming the visibilities, and a deconvolution is 
later applied to correct for the artifacts introduced by the (u, v) plane sampling. Numerous 
examples of scientific results obtained with millimeter or submillimeter interferometry can 
be found in the literature (e.g. the special edition of Astronomy and Astrophysics, volume 
468, 2007). 

To measure the source visibility on long (larger than a few hundreds meters) baselines, 
special techniques must be deployed in order to correct for the atmospheric phase errors. 
These have to be tracked on timescales of ~1 second and water vapor radiometers have 
proved to yield very efficient corrections by estimating the atmospheric contribution to the 
phase error (water content). This cophasing process is much more painful in the optical do- 
main where fringe-trackers must actually form the fringes in order to maintain the global 
wavefront as planar as possible (see e.g |Blind et al[|201 la| ). In order to obtain the best possi- 
ble images, a special attention has then to be paid to the (w, v) plane coverage that is obtained 
in such millimeter observations. Indeed, the poor sampling of the (u, v) plane has been iden- 
tified as one of the major limitations of the imaging process: the point- spread-function (PSF) 
of the synthesized instrument (the so-called dirty beam) is the Fourier Transform of the (u, v) 
coverage and may thus include strong sidelobes if the (u, v) sampling is poor, making the as- 
tronomical interpretation of the resulting images quite difficult. Although this can be partly 
corrected by the deconvolution process (CLEAN being the most widely algorithm used), 
obtaining directly a high-quality image is a central goal of the current instruments design 
and operations. The Plateau de Bure interferometer configurations have been optimized for 
that purpose ( Karastergiou et al, 2006). The positions of the 64 ALMA antennas have also 
been optimized (ALMA memo #400), the large number of antennas making it possible to 
achieve almost perfect beam shapes. In such optimization, the goal is to obtain a Gaussian 
distribution in the (u, v) plane (which means higher weights on the short baselines), which 
results in a Gaussian PSF. The behavior of a single-dish telescope with Gaussian beam taper 
is thus reproduced as accurately as possible. 
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6.3 Improving image reconstruction 

6.3.1 What is a good image reconstruction? 

In this work, we wondered whether a reconstruction is a good representation of the model. 
Although visual inspection can do a good job, it is difficult to go further. It is useful to 
compare different results from different algorithms like in the Beauty Contests ( |Lawson| 
|et al| |2004| |2006[ |Cotton et af] |2QQ8| ), but not really useful when assessing the quality of 
the image reconstruction by itself. Developing tools to compare relative positions or relative 
fluxes of the image features appears to be quite important but it is beyond the scope of this 
work. A systematic testing of image reconstruction "fidelity", such as the one mentioned 
in Appendix |B] is an essential step to build the trust in image reconstruction. The reader is 
referred to Renard et al (2011 ) for a first attempt to test systematically the quality of image 
reconstruction. 

The answer to the question of whether model fitting should be performed at the same 
time is certainly positive since this is the only way to deduce the quantitative parameters 
of the models. However, image reconstruction helps to recognize the general structure of 
the image which can be quite complex as illustrated in the supergiant and AGN cases. The 
modeler is then guided in his work by the general architecture revealed by the reconstructed 
image. 

6.3.2 Improving current facilities 

For this work, we used expected features of interferometers and instruments either in ex- 
istence or under active development. In these cases, an ideal observing sequence is often 
impossible due to technical limits: finite range of delay lines, telescope shadowing, bad 
telescope tracking near zenith, etc. Indeed, to optimally use VLTI and CHARA, we must 
identify the current limitations of instruments at VLTI and CHARA. 

VLTI. Some general remarks concerning image reconstruction at VLTI can be drawn: ( i) 
so far simple objects are the most commonly imaged: binaries, simple circumstellar environ- 
ments, or images with small field-of-view (perhaps 5x5 pixels at maximum); (ii) the main 
results have been obtained in visitor-mode with at least three configurations of 3 ATs; (Hi) 
a sampling speed of the (u, v) plane at « 25min per calibrated observation has been shown 
to be sufficient; (iv) the poor quality of the calibrated visibilities is an important limitation 
leading to poor dynamic range; (v) the impossibility to track on resolved fringes (< 10%) 
with the VLTI fringe tracking instrument FINITO is a limitation for imaging stellar surfaces; 
( vi) the observations are limited in hour angle because of the delay line stroke reducing the 
possible (u, v) coverage; (vii) finally all existing data sets suffer so far from an important 
hole in the (u, v) plane due to the current limitations in baseline configuration. In the coming 
years, the VLTI will improve imaging by prioritizing improved data quality as well as an 
increase of the number of offered interferometric baselines. Background work is on-going, 
focusing on several open questions (such as the trade-off between maximum baseline length, 
(u, v) plane density and sky-coverage). 

CHARA. The fixed locations of the CHARA telescopes mean that the array can only 
image objects of a certain size range. For near-infrared, the ideal range for imaging is be- 
tween 1.5 and 3.5 milli-arcseconds, since this allows all the 15 baselines of CHARA to 
make a useful contribution to the image reconstruction. If an object is smaller than this, then 
there is not enough resolution; if the object is bigger, then most of the short baselines are 
sampling beyond the first null and this means too much information is missing to reliably 
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invert the Fourier data. One way to improve imaging on CHARA substantially would be to 
add a telescope to the center of the arra^to introduce more small baselines and to aid the 
baseline bootstrapping. A simpler short-term solution to improving CHARA imaging is to 
use all 6 telescopes simultaneously with new combiner^] 



6.4 How to obtain high-quality images 

For future facilities, either new interferometers like MROI or new instruments on existing 
facilities, it is important to identify the most important criteria for high-quality imaging. The 
main practical limitations are related to the complexity, dynamic range, and magnitude of 
the source to be imaged. 



As pointed out in Sect. |6.2| as in the radio domain the dominant source of artifacts is 

an inadequate (u, v) sampling. The causes may come from the number of positions of the 
telescopes (see the discussions above about current limitations at VLTI and CHARA), and 
the reduced stroke of the delay lines preventing full observation of the source through the 
sky. Most specialists agree that this is the major issue by a large margin. In our simulations, 
the (u, v) coverage does not seem to be a limitation, but we have been quite generous in (w, v) 
points. Experience with the MIRC/CHARA instrument shows that snapshot coverage with 
only 4 telescopes has proven adequate for imaging objects such as stars with simple spots 
and resolved binaries if they are of the "right" size. However, more complicated objects such 
as multiple components objects, crowded stellar fields or heavily spotted stars need more 
(u, v) coverage. In addition to that, short-baseline (u, v) sampling is an essential part when 
observing extended objects (such as the evolved star case (see also |Monnier et al[ |2004)). 
Diffraction-limited imaging with a conventional single-pupil telescope, or, preferably single 
dish aperture masking, like in millimeter-wave interferometry, might be an efficient way to 
cover that issue fore relatively compact objects. 

How many telescopes are needed? Millimeter- wave interferometers have found that 6-7 
telescopes are adequate for many imaging applications and experienced observers with the 
CHARA/MIRC instrument are coming to similar conclusions. In the optical domain there 
is one crucial point that has to be integrated in this discussion: the ability to cophase the ar- 
ray. The so-called "fringe-trackers" of future generation instruments will have to cope with 
the difficulty that highly detailed images require high spatial frequency measurements. This 
translates into a strong requirement on these cophasing instruments and on the array config- 
uration: they need to be able to maintain the array phasing while measuring small visibilities. 
While this consideration was integrated in the design of NPOI (e.g. Armstrong et al[| 1998a) 
and MROI through the use of array configurations allowing for baseline bootstrapping, this 
was not the case for VLTI. The ability to re-configure telescope locations will be another 
key to making images of high spatial dynamic range, although many interesting objects such 
as spotted rotating stars and binary systems generally change in time too quickly to make 
re-configuring practical to schedule. 

Spectral dispersion has been quite successful to increase the (u, v) coverage by wave- 
length synthesis. Line maps have been obtained at NPOI for a long time (e.g. |Schmitt et al[ 
2009 ) while this has only become possible recently at VLTI with the advent of the AMBER 
high resolution spectroscopy capability Millour et al (2011). Observations in an absorption 

8 This addition is feasible since there is space on the mountain top for the telescope and space in the lab 
for the delay line. 

9 The MIRC 6-T upgrade was successfully completed in summer 2012. 
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line and in the continuum probe different atmospheric depths and thus layers at different tem- 
peratures. Such photosphere observations provide important information on the wavelength 
dependence of limb darkening and therefore possibly a strong constraint to the supergiant 
case (and all stellar- surface imaging). Another interesting application is the comparison be- 
tween a central star emitting in the continuum and its surrounding hot line-emitting disk. 
Specific algorithms could be developed to image simultaneously all channels within a line 
(e.g. taking into account the fact that a smooth transition from one channel to the next one 
is expected) although this has never been done in the millimeter wave domain where each 
channel is imaged separately, most probably because the (u, v) coverage is sufficient for such 
mappings. 

Finally, flux sensitivity and measurement accuracies appear to impose a quite stringent 
limitation. The dynamic range DN is not only a direct function of the number n of (u, v) 
points but also of the errors on each measurements 6V and 8<p as poi nted out by Bal dwin| 
and Haniff ( |2QQ2| ) who propose a rule-of-thumb estimator: DN = ^n/[(8V/V) 2 + (£0) 2 ]. 
This ideal metric relies on the independence of the measurements whose calibration should 
remain uncorrelated. Instruments should therefore be conceived to allow for precisely cali- 
brated measurements. 

The temporal fluctuations of the sources might be a severe limitation especially with in- 
terferometers having a limited number of telescope combinations. For example, the timescale 
of convection is crucial to constrain the simulations of red- supergiant stars. We probably 
need a time series of observations (at least 4-5 different epochs) about 1-2 months apart in 
the H and K bands where only small differences are predicted and at maximum 1 month 
apart in the optical where things change faster. Increasingly, we do find our prime targets 
changing on timescales of days or even hours and imaging snapshots will require as many 
simultaneous measurements as possibkpl 



7 Conclusions 

We have reviewed past and present efforts to generate aperture synthesis maps from optical 
long baseline interferometers showing the tremendous progresses obtained since 1996. We 
have presented in this paper end-to-end simulations of 4 different astrophysical cases rele- 
vant for image reconstruction with optical interferometers. The purpose of these simulations 
was to assess if reconstructed images would be faithful enough to be trusted when interpret- 
ing the results and to identify ways to increase the quality of the reconstructed images. This 
was not a trivial process and it required the contribution of a large range of competences 
from astrophysicists able to create sensible images of the objects, interferometry experts to 
mock realistic data, and specialists of image reconstruction to reach the best reconstructed 
images. We addressed four different topics of interest for optical long-baseline interferome- 
try and found the following results: 

- Detection of planet signatures in disk structures like shadows are challenged by dy- 
namic range and might be out of reach of current optical interferometry facilities. How- 
ever other larger-scale structures, with lower contrasts such as disk gravity waves, gaps, 
hot accreting bubbles around a protoplanet etc. might be revealed. Future efforts in that 
direction should address both instrumental issues (precision and accuracy of measure- 
ments, data reduction) and image-reconstruction dynamic range limits. 

10 A warning to researchers that believe building up many 3-telescope observations is the best approach: it 
takes 20x longer to build up full closure phase coverage one triplet at a time for a 6-telescope interferometer 
compared to schemes which measure all baselines!. 
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- Interferometric imaging allows sizes of convection cells and the contrast of the surface 
intensity to be measured and allows us to better understand the photosphere and dynam- 
ics of evolved massive stars. These results are crucial for understanding their unknown 
mass loss mechanism that contributes heavily to the chemical enrichment of the Galaxy. 
The complex structure requires a rich (u, v) plane sampling and therefore the combina- 
tion of the highest possible number of telescopes. 

- The capabilities offered by current interferometers make it possible to reveal asymme- 
tries in the molecular layers around evolved stars. Further imaging of molecular layers 
appearing in spectral lines is the path to understand these dust and molecule factories of 
the Galaxy. 

- Our simulations show that unraveling the clumpy structure of the dusty torus obscuration 
region in AGN is within reach and will certainly be a necessary step to progress in the 
understanding of the AGN unified scheme. 

In general, we found that image-reconstruction techniques have been successful at retrieving 
global morphologies of all our objects. Their current limitations lie (i) in their ability to 
match exact intensities and (ii) to detect highly contrasted features. 

The (u, v) coverage is of crucial and primary importance and therefore we look forward 
to results from the 6-telescope mode of CHARA/MIRC, the implementation of at least 6 (up 
to 8) telescopes at the VLTI with a fully operational set of delay lines, and the construction 
of MROI with full capacity. We also recall the importance of imaging instruments (Matisse, 
Gravity, VSI at VLTI, MIRC and MIRC6T at CHARA). The spectral resolution is an asset 
whose full exploitation has only started to be fully exploited. Of course, like in millimeter 
astronomy, the future of interferometry has to be prepared and will probably aim at a facility 
similar to the Very Large Array (VIA) or ALMA. 

Our study also shows that, in the context of image reconstruction, the tools have reached 
a certain level of maturity to produce on-the-fly results and with a reliability and repro- 
ducibility that allows productive results and discussion. Optimization of the "theory — » 
models — > simulated observations — > reconstructed images" cycle is never finished, but a 
first level of quantitative analysis shows that optical interferometry is capable of produc- 
ing meaningful images. We have now tools to improve the comprehension of astrophysical 
phenomena and to reach ambitious research topics. 

We believe that imaging interferometry will provide to the broad astronomical commu- 
nity a new window to the universe if we are able to gather resources (expertise, manpower 
and relatively limited investments). 
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A Noise model for ASPRO simulations 

This Appendix is aimed at describing the noise model used in ASPRO for creating the simulations used in 
this paper (see Sect. |5.2) . 

The noise model is based on a general scheme valid for spatially filtered recombiners where the detection 
of fringes is made on a detector with "pixels". This scheme is valid for image-plane recombination, with 
fringes covering a surface of a pixel camera, and for pupil recombination where fringes are obtained on a few 
pixels detector by scanning in optical path difference. 

The flux A^t from the object of magnitude ra& in a given bandwidth AA of a photometric band b is 
collected by iV te i telescopes of diameter D, transmitted with some instrumental transmission 7^, and injected 
with some Strehl factor s due to incomplete correction of wavefront aberrations due to seeing into a spatial 
filter like an optical single mode fiber for example, during a time t mt . Thus: 

N^ = rjF lO-° Amb T imt N tel s ^ AA, (5) 

where F$ is the zero-magnitude flux in band b, expressed in phs -1 m _2 /mi _1 transmitted through the atmo- 
sphere with an absorption 77. 

This flux is divided between the photometric flux and interferometric flux with a branching value bi, 
where bi equals 1 for recombiners which do not have simultaneous flux monitoring. 

The N te \ photometric fluxes N p = (l-bi) N tot /N te \ are distributed on N p - lx pixels. The interferometric flux 
Ni = bi Mot /Mel consists in Nf = N te \ (N te \-l)/2 fringes that cover pixels. There are thus N pp f = N^JNf 
pixels per fringe. 

These fringes code the intrinsic visibility V(u, v, A) degraded by the interferometer instrumental contrast 
and the atmosphere (through the jitter associated with the temporal coherence of the seeing). V(u, v, A) and the 
derived interferometric observables are thus affected by the sum of the variance of the flux used to code the 
corresponding fringe in the interferometric flux and of the associated two photometric fluxes. For example, 
since the squared visibility estimator of a correlated flux F l c J measured alongside with photometries F; and 
Fj is V 2 = \ < \F l c^ > I < F t Fj >, the associated variance is 

o- 2 (V 2 ) = NiV(u, v, A) + N pp fcr 2 et + 2(N~ p + A^det) , (6) 

where cr^ is the readout noise of the detector. 

The noise model used in ASPRO takes also into account the possibility of increasing the integration time 
to keep observations in a photon-dominated regime, when a fringe tracker is present. 

Finally, no detailed calibration error was computed. We took instead an additional visibility and closure 
phase threshold error set to 0.002 in visibility, and 0.1 degree in closure phase. 



B Computing fidelities 

Testing image reconstruction software is out of scope of this paper. In this Appendix, we follow the approach 
described in ALMA Memo 398 (F. Gueth private communication) to evaluate the quality of the reconstructed 
images presented in this paper. One of the possible methods to compare original (convolved to the interfer- 
ometer resolution) and reconstructed image is to compute the fidelity of the image. This can be done either 
in the direct image plane or in the spatial frequency (u, v) plane. Such a pixel to pixel comparison requires 
subpixel alignment to limit the effect of sharp transitions. 
In the image plane this fidelity can be expressed as 

T(x y) = abs(MoJg/(x,y)) ^ 
max(Diff(x, y), Threshold) ' 



11 Available at http://www.jmmc.fr/aspro 



Available at http://www.maumae.net/yorick/doc/index.php 
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Table 6 For a given fraction of the brightest pixels (col 1) the corresponding number of pixels, median and 
average fidelity is given for each object. 





Supergiant 




Evolved Star 






AGN 




(100 x 100 pixels) 


(200 x 200 pixels) 


(200 x 200 pixels) 


Fraction of 


Npixel 


Med 


Avg 


Npixel 


Med 


Avg 


Npixel 


Med Avg 


brightest pixels 


















0.3% 


30 


2.9 


3.2 


119 


13.7 


30.7 


122 


1.4 2.9 


1% 


100 


3.6 


4.8 


398 


4.4 


14.0 


408 


1.4 4.5 


10% 


1000 


5.7 


11.5 


3980 


0.0 


1.4 


4080 


1.5 4.5 



where Model(x,y) describes the object "true" brightness distribution and 

Diff(x, y) = Model(x, y) - Reconstructed(x - Ax, y - Ay) (8) 

describes the difference between the model and the reconstructed image shifted by the offset (Ax, Ay) to have 
the images centered. For a proper comparison images are normalized to the total intensity contain in the 
image. Threshold is defined here as 0.7 rms[Diff(x, y)] which provides an estimation of the threshold noise of 
this difference. For example a pixel fidelity of 100 corresponds to a difference of 1% between the model and 
reconstructed image. We have computed such a fidelity for three of our objects: the supergiant, the evolved 
star and the active galactic nuclei. Figures Qo| and offer different ways to visualize fidelity. The first row 
shows the cumulated average fidelity over the image for the three different objects. The second row displays 
the original model image filtered by fidelity level, i.e. only the pixels with fidelity above a certain level are 
displayed. Finally Table[6]shows the average fidelity of the pixels whose intensity is above a certain fraction 
of the total image intensity. 
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Fig. 10 Top row: cumulated average fidelity for the supergiant and evolved star cases. Blue dashed, green 
dashed-dotted, red dotted vertical lines: number of pixels with intensity greater than a given fraction of the 
total image intensity, respectively: 4.10 -4 , 1.10 -4 , 4.10 -5 for the supergiant (left), 5.10 -3 , 1.10 -3 , 5.10 -4 for 
the evolved star (right). Horizontal lines share the same color code: median fidelity for pixels having intensity 
greater than a given fraction of image total intensity. Central row: original model image filtered with a 
fidelity threshold (supergiant: 10, evolved star: 10). Bottom row: fidelity profile along the solid cut line in 
central row. Orange solid line: fidelity; green dashed-dotted line: model image; dashed blue: reconstructed 
image. Units are given in fidelity, image profiles have been scaled to fit the figure. 
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Fig. 11 Top: cumulated average fidelity for the AGN case. Blue dashed, green dashed-dotted, red dotted 
vertical lines: number of pixels with intensity greater than a given fraction of the total image intensity i.e. 
1.10 -3 , 5.10 -4 , 1.10 -4 . Horizontal lines share the same color code: median fidelity for pixels having intensity 
greater than a given fraction of image total intensity. Central: original model image filtered with a fidelity 
threshold (AGN: 3). Bottom: fidelity profile along the solid cut line in central row. Orange solid line: fidelity; 
green dashed-dotted line: model image; dashed blue: reconstructed image. Units are given in fidelity, image 
profiles have been scaled to fit the figure. 



